#include <math.h>

//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
// Function:  InvNorm
//
// Purpose:   Inverse of normal distribution, returning x(y) where
//                 y(x) = Normal(x)
//                           x
//                      = INTEGRAL exp(-u*u/2) du / sqrt(2*pi)
//                         -inf
//                      = ( 1 + erf(x/sqrt(2)) ) / 2
//
//            Rational approximations claimed to be correct to 1.15e-9
//
// History:   JS        13 Sep 2003
//            From Peter J Acklam www//home.online.no/~pjacklam/notes/invnorm
//-----------------------------------------------------------------------------
//
double inv_norm(   //   O Inverse normal (= number of standard deviations)
				double  y)        // I   Argument       (= cumulative probability)
{
	double  r, s;
	if( y < 0.02425 )
	{
		s = sqrt(-2.0 * log(y));
		return  + (2.938163982698783  +
			(4.374664141464968  -
			(2.549732539343734  +
			(2.400758277161838  +
			(0.3223964580411365 +
			0.007784894002430293 * s) * s) * s) * s) * s)
			/ (1.0                +
			(3.754408661907416  +
			(2.445134137142996  +
			(0.3224671290700398 +
			0.007784695709041462 * s) * s) * s) * s);
	}
	else if( y > 0.97575 )
	{
		s = sqrt(-2.0 * log(1.0 - y));
		return  - (2.938163982698783  +
			(4.374664141464968  -
			(2.549732539343734  +
			(2.400758277161838  +
			(0.3223964580411365 +
			0.007784894002430293 * s) * s) * s) * s) * s)
			/ (1.0                +
			(3.754408661907416  +
			(2.445134137142996  +
			(0.3224671290700398 +
			0.007784695709041462 * s) * s) * s) * s);
	}
	else
	{
		s = y - 0.5;
		r = s * s;
		return  s * (2.506628277459239 -
			(30.66479806614716 -
			(138.3577518672690 -
			(275.9285104469687 -
			(220.9460984245205 -
			39.69683028665376 * r) * r) * r) * r) * r)
			/ (1.0 -
			(13.28068155288572 -
			(66.80131188771972 -
			(155.6989798598866 -
			(161.5858368580409 -
			54.47609879822406 * r) * r) * r) * r) * r);
	}
}